MICROCOPY  RESOLUTION  TEST  CHA5?T 

NAriONAl  BURtAIJ  OF  SMNDAffOS  I9b3-^i 

\ 


NPS-69-78-012 


NAVAL 


POSTGRADUATE  SCHOOL 

Monterey,  California 


D D c 


pjmrsEiiD.r'j? 


JUL  18  1978 


fEISEETTIL 


o 


A CCMUTER  SUBROUTINE  FOR  STRESS  ANALYSIS  OF 
ROTATING,  HEATED  DISKS 

by 

John  E.  Brock 
Robert  E.  Brown 
May  1978 


Approved  for  public  release;  distribution  unlimited. 


Prepared  for: 

Chief  of  Naval  Research 
Arlington,  Virginia  22217 


78  07  10  019 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


Rear  Admiral  Tyler  Dedman 
Superintendent 


A COMPUTER  SUBROUTINE  FOR  STOESS  ANALYSIS 
OF  ROTATING,  HEATED  DISKS 


This  report  gives  listing  and  instnx:tions  for  using  a 
digital  conpiter  subroutine  for  finding  stress  distri- 
bution in  a thin  rotating  disk  with  nonuniform  heating; 
the  problem  is  axisymnetric.  An  iterative  method  is  used. 
Theoretical  background  is  given. 


J.  R.  Bursting 
Provost 


Approved  by: 


Allen  E.  Fuhs,  Chairman  W.  M.  Tolies 

Mechanical  Engineering  Department  Dean  of  Research,  Acting 


UNCLASSIFIED. 


SteCIHTV  CLASSIFICATIOM  OF  THIS  FACE  (Whtn  Dtm  gin»f»«p 


hps- 


REPORT  DOCUMENTATION  PAGE 


READ  mSTRUCnC  Nt 
BEFORE  COMPLETWC  TORM 
3.  rccificnt's  cataloo  MUMatn 


FQIIT.liUMBU- 


NPS-  69-78-jrL2 


2.  GOVT  accession  NO. 


S.  TYPE  OF  REPONT  • PCMOD  COVCNCO 


^ CCMPUTER  SUBROUTINE  FOR  STRESS  ANALYSIS  OF 
rotating,  H^TED  DISKS-  ^ 


».  PERFORMING  ORG.  REPORT  NUMBER 


K 


t.  CONTRACT  OR  GRANT  NUMRCRT*) 


John  E. /Brock 
Robert  E. /Brown 


’ w¥:"ffr5a‘  t!sai°!SBc) 

Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 


10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  • WORK  UNIT  NUMBERS 

611 52N,  RROOO-01-01 
N0001478WR80023 


It.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Naval  Postgraduate  School  ^ 
Monterey,  California  93940 


REPORT  DATE 

' May  iR78  J 


AGES 


^4.  MONITORING  AGENCY  NAME  A ADDRESS/// from  Controlling  OUic») 

Chief  of  Naval  Research 
Arlington,  Virginia  22217 


IS.  SECURITY  CLASS.  (Mt^raperi; 


/ 


1S«.  OECLASMlCAXi^ 
SCHEDULE 


■ RAOING 


16.  Distribution  statement  (oi  thim  Rtport) 


Approved  for  public  release;  distribution  unlimited. 


17.  distribution  statement  (ot  Ihm  mbttrmel  tntmd  In  Block  30.  if  dlllotcni  Item  Report) 


■I 


IB.  SUPPLEMENTARY  NOTES 


It.  KEY  WORDS  (Cenllnue  on  rerereo  elde  It  noceoeerr  mtd  Idonllty  ky  block  nuaibarj 

Stress  analysis 

Rotating  disks,  Heated  Disks 

Axisyiimetric  Elasticity,  Axisymmetric  Disks,  Elastic  Disks 


'^^mAcT7ConUnuo'o!rtoretoo'oldirillnocoeell!r‘mdldoniiir'^TioeknimSer)'~~~~~~^~~~~~~~'~^~~~~~~~~ 

This  report  gives  listing  and  instructions  for  using  a digital  conqjuter 


20.  A 


subroutine  for  finding  stress  distribution  in  a thin  rotating  disk  with 
nonuniform  heating;  the  problem  is  axisymmetric.  An  iterative  method  is 
used.  Theoretical  background  is  given. 


1 

I 


\ 


DD  , 


FORM 
JAN  73 


1473 


EDITION  OF  I NOV  61  IS  OBSOLETE 
S/N  0I02*014>6601  I 


.^OBSOLETE  ^ , ry  UN(S4pilFIEf) 

d \AeAriTV  CL'TOIFTCATioN  of  this  page 


fWhwi  Dele  Bnlorod) 


TABLE  OF  CONTENTS 


Introduction — 1 

Fundamental  Analysis - 2 

Power  Law  of  Thickness  Variation 3 

Exponential  Law  of  Thickness  Variation-- 4 

General  Law  of  Thickness  Variation — 5 

Computer  Implementation 8 

References 9 

Appendix  A - — - 10 

Appendix  B-- 13 

Appendix  C 16 

Initial  Distribution  List 20 


■Ami  I 


. — ^ 


- J-  - 


A Ccxnputer  Subroutine  for  Stress 
Analysis  of  Rotating,  Heated  Disks 

by 

John  E.  Brock  emd  Ptobert  E.  Brovm 

Introduction 

Although,  as  is  indicated  by  the  title  hereof,  the  principal 
purpose  of  this  monograph  is  to  present  tested  and  proved  digital  cran- 
puter  software  for  the  analysis  of  stress  in  a spinning  axisyirrvetric 
disk  having  a radially  variable  thermal  strain  field,  the  opportunity 
is  cilso  taken  of  developing  the  theory  and  oresenting  seme  analytic 
solutions . 

The  method  developed  herein  for  computer  analysis  of  disks 
having  a general  law  of  thickness  variation  was  suggested  by  the  al- 
gorithm contained  in  reference  2 and  it  appears  to  have  advantages 
over  such  procedures  as  that  of  M.  Donath,  reference  3>  which  Pias  been 
widely  circulated  in  a book  by  S.  Timoshenko,  reference  5. 

In  what  follows  we  tanediately  obtain  a second  order  linear 
differential  equation  with  dependent  variable  u,  the  radial  defomatlon, 
and  r,  the  radius.  Analjrtic  treatment  is  given  for  two  particular  laws 
of  thickness  variation.  For  the  general  case  of  thickness  variation, 
the  equation  is  recast  as  a second  order  linear  differential  equation 
in  which  the  dependent  variable  is  radial  stress,  o^.  However,  for 
numerical  treatment  an  alternate  form  is  more  useful  and  direct  and 
this  forms  the  basis  of  the  digital  computer  software  which  is  given 


n 


I 


f 


and  illustrated  in  the  appendices  hereto. 

Fundamental  Analysis 

We  presume  that  the  disk  Is  thin  enough  and  that  the  thick- 
ness varies  slowly  enouggh  with  respect  to  radius  that  we  are  Just- 
ified in  neglecting  all  stress  components  excepting  only  the  radial 
stress  Op  and  the  circumferential  stress  Og.  Material  properties 
E,  Young's  modulus,  and  v,  Poisson's  ratio,  are  presumed  to  be  in- 
deed constant.  The  thermal  strain  field,  oT,  and  the  density  y 
may  be  specified  functions  of  radius. 

Two  types  of  problem  are  considered: 

1.  Annular  disk,  0 < a S r = b,  with  oAa.)  and  a^(b) 

r r 

being  specified. 

2.  Solid  disk,  0 a r = b,  with  a^Cb)  being  specified. 

We  also  use  the  symbols  t = t(r)  for  disk  thickness  and  u for 
angular  velocity.  Other  simplifying  notation  will  be  introduced 
later  on. 

Consideration  of  radial  equilibrium  leads  without  difficulty 


to  the  equation 


1 d(rta  ) 
t - dr  ■ 


= Og  - 


(1) 


The  thermoelastic  constitutive  equations  are 

Eg-  » Oq  - va  + BaT;  Eg^  = c - vo-  + BaT  (2a, b) 

u o r r r “ 

where  the  strain  con^xsnents  are 

Gg  = u/r;  Gj.  = du/dr  (3a,b) 

Strain  corpatabillty  leads  to  the  equation 

r ^Eu/r)  * -(l+v)(ag-Op)  (*4) 

These  equations  njay  easily  be  combined  into  the  differential  equation 
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u"  + u'/r  - u/r*  + (tVt)[u'  + vu/r  - (l+v)aT] 

■ (l+v)(oT)'  - (l-v*)Yto^r/E  (5) 

where  primes  denote  differentiation  with  respect  to  r. 

Two  particular  laws  of  thickness  variation  permit  slmole 
analytic  treatment. 

Power  Law  of  Thickness  Variation 
If 

t - to(r/ro)”  (6) 

so  that 

(t'/t)  = n/r  (7) 

then  equation  5 becomes 

u"  + (l+vn)u'/r  + (vn-l)u/r^  = 8*  + nS/r  - kr  (8) 

where 

8 =*  (l+v)aT,  k = (1-v^)yw^/E  (9,10) 

Equation  8 may  be  rewritten  as 

„)]  = ^(2+n-m)/2  ^ 

where 

m = ±/(n^-A  n+'4)  » ±/C(n-2)^+i)(l-  )n]  (12) 

and  either  the  positive  or  the  negative  sign  may  be  used  Equation  11 
may  be  proved  simply  by  performing  the  indicated  operations  and  com- 
paring with  equation  8. 

The  quantity  on  the  right  in  equation  11  is  well  defined  so 
that  the  solution  of  the  differential  equation  may  be  obtained  slnoly 
by  integration,  multiplication  by  r"'“^,  and  another  integration.  Two 
constants  of  integration  are  introduced.  For  the  solid  disk  (case  2) 

- 3 - 


r 


1 

i 

I 

I 

I 


I 

1 

1 


1 

i 


u(0)  = 0 gives  one  condition  and  the  second  ccxnes  from  satisfying 
the  given  value  of  a^Cb).  For  the  annular  disk  (case  1)  satisfying 
the  given  conditions  c^Ca)  and  a^(h)  permits  evaluating  the  con- 
stants. 

An  exc  ':.e  with  n = - 0.42  ,3  given  in  Appendix  C.  ^fote 
that  the  case  n “ 0 corresponds  to  a disk  of  uniform  thickness. 

Exponential  Law  of  Thickness  Variation 
If 

t = to  exp(-mr^)  (I3) 

where  m is  a constant  of  appropriate  dimensionality,  then 

(tVt)  ® -2rm  (li4) 

and  equation  5 becomes 

u"  + (1/r  - 2rm)u’  - (1/r^  + 2vm)  = -2m8  S'  - kr  (15) 
If  additionally  we  assume  that  8 ' = 0 (which  makes  the  themal 
strain  field  constant  — a triviality)  the  solution  is  simply 

u = r(3+k/2m)/(l+v)  (I6) 

In  this  case  we  can  easily  find 

cfp  = O9  = YuV2m  (17) 

which  is  independent  of  r.  Thus,  if  an  allowable  normal  stress 
is  specified  and  if  blade  or  bucket  loading  at  r = b is  w (oounds, 
say)  per  unit  circumference,  then  a disk  having  thickness 

t = (w/a^)  exp[(b^-  r^)(Y&)V2a^)]  (I8) 

will  be  such  that  = Cg  h a^.  If  the  failure  criterion  is  the 
maximum  shearing  stress  criterion  (Tresca's  condition),  it  is  clear 
that  this  disk  is  optimal  in  the  sense  of  having  least  volume  and 
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thus  having  least  weight.  If  the  failure  criterion  were  that  of 


von  Mises,  a slightly  lighter  disk  would  suffice. 

General  Law  of  Thickness  Variation 
Equations  3a  and  2a  give 

Eu  = rBaT  + r(ag-vOp)  (19) 

au^  by  differentiation 

Eu'  =>  EhT  + rE(aT)'  + (Og-vOp)  + r((jg’-v0p’)  (20) 

Ptoti  3b  and  2b  we  also  have 

Eu'  = BaT  + - vCfl  (21) 

r 0 

Subtracting  21  from  20  and  rearranging  gives 

E(aT)'  + (l+v)(Og-o^)/r  + Og'  - * 0 (22) 

Ev^uation  1 may  be  rewritten  as 

Og-Op  3 ro^'  + rvOp  + yw^r^  (23) 

where  we  have  written 

V » t'/t  (2^) 

for  convenience.  Differentiating  23  we  get 

o-'  » 2o  ' + ro  " + vo^  + rva'  + rv'a^  + 2Y&>^r  (25) 

0 r r r r r 

and  substituting  23  and  25  into  22  gives 


ra  " + (3+rv)o  ' + [(2+v)v+rv' ]o_  + (3+v)Yu*r  + E(aT)’  = 9 (26) 
r r r 

This  is  a single  differential  equation  with  dependent  var- 
iable 0 and  can  be  dealt  with  by  standard  numerical  methods.  The 
r 

conditions  for  evaluating  the  constants  of  integration  have  been 
mentioned  earlier. 

However,  the  preceding  procedure  is  not  particularly  satis- 
factory. For  one  thing,  the  solid  disk  case  for  which  r can  become 
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zero  encounters  numerical  difficulties  unless  special  treatment  is 
employed  to  evade  them.  More  importantly,  however,  if  T is  given 
by  a graph  or  a numerical  table,  determination  of  (oT)'  may  Involve 
rumerlcal  difficulties  which,  ultimately,  are  due  to  our  having  per- 
formed the  differentiation  to  arrive  at  equation  20.  Accordingly, 
an  alternate  procedure  which  adheres  more  closely  to  the  fundamental 
mechanics  of  the  problem  is  described  in  what  follows. 

We  consider  the  annular  disk  first  and  we  represent  the  un- 
known stress  difference  in  the  form 

Og-Op  = (A+Bn)r  (27) 

where  A and  B are  unknown  constants  and  n is  an  unknown  function 
normalized  so  that 

n(a)  = 0,  n(b)  » 1 (28a, b) 


Initially  we  make  an  assumption  for  n taking  a linear  var- 
iation in  the  absence  of  better  information.  By  the  use  of  some  of 
the  preceding  equations  we  will  be  able  to  construct  am  improved 
form  for  n smd  will  iterate  until  there  is  satisfactory  convergence. 
Let 

z = Eu/r,  w » BaT,  B ® (29,30,31) 

noting  that  w aind  S have  different  meamlngs  here  than  when  they 

were  used  earlier.  We  cam  recaist  equation  1 as 

d(to„)/dr  ■ (aQ-a^)t/r  - yt^^rt  = t(A+Bn-Yaj^r)  (32) 

r y r 

Before  performing  the  indicated  integration  we  introduce  two 


convenient  notatloral  devices,  viz. 

p...  « f ...  dr;  *...  = C...X 

Jo 


(33, 3^*) 


Thus,  from  equation  32  we  have 
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ta  = + Apt  + Bpnt  - pBrt 

r*  r*  « 

From  equation  4 and  equation  27  we  have 
dz/dr  » -(l+v)(A+Bn) 


so  that 


Equation  2a  is 


so  that 


z = (z)  - (l+v)CA(r-a)  + Bpn] 

A 


z = w + (o.-a)  + (l-v)o^  = w + (A+Bri)r  + (l-v)a  (38) 
y r r r 


(z)  = (w)  + Aa  + (l-v)((j  ) (39) 

cl  ^ r 21 

(z)^  = + (A+B)b  + (l-v)(a^)^  (40) 

Evaluating  38  at  r » b and  using  39  and  40  gives  one  equation 

involving  the  unknowns  A and  B.  Evaluating  35  at  r = b gives  a 

second.  These  equations  can  be  airanged  as 

*(pt)  *(prit)  ■ A'  _ ■»(p8rt)  + (tap)^^  - 

(2+v)(b-a)  b+(l+v)»(pn)  _B  .(w)^-(w)^+(l-v)[(a^)^-(ap)^' 

so  that  one  can  easily  solve  for  A and  B.  Then  is  obtained  from 

35,  z is  obtained  from  37  and  39,  and  (sJq-^j.)  and  Og  are  obtained 

from  38.  Then  a new  function  n is  calculated  from 

n » [(CTg-Cp)/r-A]/B  (42) 

Using  the  new  n the  process  is  iterated,  convergence  being 

monitored  by  examination  of  the  sequence  of  values  of  A and  B that 

are  calculated.  When  convergence  is  satisfactory,  the  desired 

functions  a and  Oa  are  at  hand, 
r y 

The  situation  is  simpler  for  the  solid  disk,  case  2.  (cr^)^ 

is  not  given  but  conditions  of  continuity  require  that  (ag-Op)/r 
vardsh  at  r = 0.  Thus  A is  zero  and  B can  be  obtained  flxim 


L. 
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The  ancillary  subroutines  manipulate  these  vectors.  All  of  these 
are  obvious  except  possibly  IITTV  which  performs  an  intec^tion  by 
use  of  Milne's  formulas,  cf.  reference 

In  the  subroutine  P?0DISK  there  is  a slight  departure  from 
the  theory  as  given  herein.  As  a first  step,  all  quantities  and 
functions  were  "dedimenslonalized"  but  otherwise  the  method  is 
Just  as  described  above.  Somewhat  finer  details  of  what  was 
actually  done  may  be  found  in  reference  1. 

Appendix  C contains  some  examples  and  remarks  concerning 

them. 
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Appendix  A 

Listing  of  subroutine  RODISK 

(The  listing  on  this  page,  page  10,  is  of  the  comnents  ’^ich 

provide  instructions  for  the  use  of  RODISK.  Cormands  appear 

on  the  following  two  pages . ) 

SUBROUTINE  RODISK.  JOHN  £.  9RGCK.,  1 MAY  1978 
THIS  IS  A SUBROUTINE  FOR  DETERMINING  RADIAL  AND  CIRCUM- 
FERENTIAL STRESSES  IN  AN  AXISYMMETRIC  THIN  ELASTIC  DISK, 
ROTATING  A-^  ANGULAR  VELOCITY  CMEGA  (RADIANS  PER  SECOND) 

ABOUT  ITS  AXIS  OF  SYMMETRY  AND  HAVING  AN  AXISYMMETRIC 
DISTRIBUTION  CF  THERMAL  STRAIN.  TWO  TYPES  OF  PROBLEM  MAY 
Be  TREATED: 

TYPE  1.  ANNULAR  DISK  0*=  INSIDE  RADIUS  ARAD  AND 
OUTSIDE  RAD'^US  BRAD.  THE  RADIAL  STRESS  IS  SRA  AT  THE 
INNER  RADIUS  AND  SR6  AT  THE  OUTER  RADIUS.  THE  INSIDE 
RADIUS  MUST  PE  GREATER  THAN  ZERO. 

TYPE  2.  SOLID  DISK  OF  OUTSIDE  RADIUS  BRAD  AND  WITH 
RACIAL  STRESS  SRB  AT  THE  OUTSIDE  RADIUS. 

THE  USER  MUST  PROVIDE  A MAIN  PROGRAM  WHICH  CALLS  SUBROUTINE 
RODISK  AFTER  IT  HAS  SUPPLIED  THE  FOLLOWING  INFORMATION. 

(1)  N,  INTEGER.  (N-l)  IS  THE  NUMBER  OF  EQUAL  SUBDIVISIONS 
INTO  WHICH  the  ANNULAR  RADIUS  (BRAD  MINUS  ARAD)  IS 
DIVIDED  FOR  COMPUTATIONAL  PURPOSES.  t^e  PRESENT 
DIMENSIONING  CAN  ACCOMMODATE  N NOT  LARGER  THAN  101. 

(2)  BRAD 

(3)  ARAD  (NOT  NECESSARY  FOR  PROBLEMS  CF  TYPE  2.) 

(4)  SRB 

(5)  SRA  (NOT  NECESSARY  FOR  PROBLEMS  OF  TYPE  2.) 

(6)  TEEBEE,  DISK  ‘I’HICKNESS  AT  OUTSIDE  RADIUS 

(7)  POrS,  POISSON'S  RATIO. 

(8)  KP(1)  = 1,2.  INTEGER  TO  DENOTE  PROBLEM  CF  TYPE  1,2. 

(9)  KP(2),  integer  TO  PROVIDE  FOP  SKIPPING  WHILE  PRINTING 
OUTPUT.  FOR  EXAMPLE,  IF  N=101  AND  KP(2)=5,  ONLY 
EVERY  FIFTH  SET  OF  VALUES  WILL  BE  PRINTED:  1ST,  6TH , 
...,  96TH,  lOlST. 

(10)  KP(3),  INTEGER  SPECIFYING  THE  NUMBER  OF  ITERATIONS 
TO  BE  PERFORMED.  USUALLY  KP(3)=10  IS  SUFFICIENT  FOR 
ENGINEERING  ACCURACY. 

(11)  KP(4).  IF  KP(4)=0  ONLY  FINAL  ANSWERS  WILL  EE  PRINTED. 
IF  KP(4)=1  a SEQUENCE  OF  ITERANT  VALUES  WILL  BE 
PRINTED,  INDICATING  DEGREE  OF  CONVERGENCE. 

(12)  VECTORS  X(I,J),  1=1,2, 3;  J= 1 , 2 , . . . ,N. 

VECTOR  X(1,J)  CONTAINS  VALUES  OF  THE  RATIO  (LOCAL 
THICKNESS  Tip  OISK)/(TEE0EE)  COMPUTED  AT  EQUALLY  SPACED 
RADII  STARTING  AT  THE  INSIDE  AND  ENDING  AT  THE  OUTSIDE, 

VECTOR  X(2,J)  CONTAINS  VALUES  OF  ( GAMMA) (OMEGA-SQUARED ) 
TIMES  ( BP.AO-S&ARED)  DIVIDED  BY  (SRB).  FOR  MOST  PROBLEMS 
GAMMA  DOES  NOT  VARY  WITH  RADIUS  AND  THIS  QUANTITY  IS  A 
CON STANT. 

VECTOR  X(3,J)  CONTAINS  VALUES  OF  ( E) ( ALPHA ) (TEE )/ (SRB ) 
WHERE  (E)  IS  YOUNG'S  MODULUS,  (ALPHA)  IS  THE  COEFFICIENT  OF 
LINEAR  THERMAL  EXPANSION,  AND  (TEE)  IS  TEMPERATURE  CHANGE, 
THE  main  PROGRAM  MUST  CONTAIN  TH=  STATEMENTS: 

IMPLICIT  REAL*8  (A-H,0-Z) 

INTEGER  KP(4) 

COMMON  X,N 

COMMON  /ONE/ARAO, BR AD ,SRA , SR B .TEEBEE , PO IS ,KP 

FOLLOWING  SUBROUTINE  RODISK  THERE  ARE  SEVERAL 
ANCILLARY  SUBROUTINES  WHICH  PERFORM  VARIOUS  OPERATIONS 
ON  THE  VECTORS  X(I,J).  THE  FUNCTION  OF  EACH  IS  CBVIOUS 
FROM  THE  LISTING.  THEY  MAY  BE  USED  IN  THE  USER'S  MAIN 
PROGRAM.  TWO  CF  THESE  ANCILLARY  SUBROUTINES  WHICH  ARE 
AVAILABLE  IN  THIS  PACKAGE  8LT  WHICH  ARE  NOT  CALLED  BY 
SUBROUTINE  PODISK  ARE  OUPV  WHICH  DUPLICATES  A VECTOR  AND 
PRIV  WHICH  PRINTS  a VECTOR. 
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SUB'^OUTIN?  ROOISK 
IMPLICIT  '^£AL*3(A-H,0-Z» 

REALMS  X(20,101  ) 

INTEGER  KP(4) 

COMMON  /6nE /ARADf  brad  ,SRA,SRBfTEE0EEfFOIStKP 

aNE=I  .D+0 

P0I5=3.0-1 

ZER3=C.0+0 

IF(KP(1>.EQ.2)  ARAD=ZER0 

RHC=ARA0/BRA0 

ENM=N-1  . . 

WRITE(6,2)  KP(l) 

2 FORMAT!//, lOX, 'RODISK  PROBLEM  OF  TYPE  *,11, 
DO  5 1=1, N 
5IM=I-1 
Y=EIM/ENM 

X(A,I  )=kHO+(ONE-RHO)*Y 
X (5,n=Y 
5 X (6,1  )=Y 
I T£R=1 

ETA=X(1,1)/X(1,N) 

IF(XP(1),E0.2)  GO  TO  100 

THE  PROBLEM  IS  OF  TYPE  OME:  ANNULAR  DISK 
SRAT=SRA/SRB 

Cl=(2.D+0+P0I Sj^IONF-RHO) 

CALL  INTV(1,7» 

C2=X(7,N) 

C5=X(3,l)-X(3,N)-(ONE-PCIS»*(ONE-SRAT ) 

CALL  MULV(1,2,8) 

CALL  MULV(8,4,9) 

CALL  INTV(9,10) 

C6=X(1C,N)+(0NE-ETA*SRAT)/(0NE-RH0) 

20  CALL  INTV(6,ll) 

C3=0NE+(0NE-RH0)=(0ME+P0!S»«X(11,N) 

CALL  MULV(1,6,12> 

CALL  INTV(12,13) 

CA=X(  13, N» 

0=C1*CA-C2-C3 

A = (C5^C4-C6=«'C3)  /O 

B=(C1*C6-C2*C5)/0 

IF(KP(4).EQ.ll  ViRITE(6,7)  ITER,  A,  3 
7 FGRMAT(5X,I 10,1P2E20.51 
CALL  MULS(7,14,A) 

CALL  MULS(13,15,B) 

CALL  A0DV( 14, 15,15) 
call  SUBV( 15, 13,16) 

S=CNE-RHO 

CALL  MULS( 16, 16,S) 

S = ETA=»SRAT 

CALL  A0DS(16,16,S) 

CALL  DIVV( 16, I, 20) 

ZJ=A*RH0+X( 3, 1 )+ (ONE-PCI S )*SR AT 
CALL  MULS( 11,11 ,3) 

CALL  MULS(5,16,A) 

CALL  ADDV(H,16,16) 

S=-(ONE-RHO)*(OKc+PDI S) 

CALL  MULS( 16, 16,S) 

CALL  A00S(16,19  ,Z0) 

S=POIS-aNE 

CALL  MULS(20, 18,S) 

CALL  ADOV( 18,19,18) 

CALL  SUBV( 18,3, 19) 

ITER=ITER+1 

IF(ITER.GT.KP(3))  GO  TO  200 
CALL  DIVV(  19,4,  I T) 

CALL  S'JBS(17,17  *) 

CALL  CIVS(17,6,bi 
GO  TO  20 
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100  ci*pcis-on: 

C2=«ETA 

C5=X(3tl)-X 13,N)+POIS-rNP 
CALL  PULV(l,2»7) 

CALL  MULV(7,4,8) 

CALL  INTV(a»9) 

C6=3Nc+X(9,N) 

105  CALL  INTV(6,10) 

C3=ONE+(ON£+POIS)*X( 10»N) 

CALL  PULV(1,6,11) 
call  INTV(11,12) 

C4»X( 12, N) 

0*C1«CA-C2=>'C3 
SPAT*  (C5*C9-C6*C.3)/0 
B*(C1<‘C6-C2*C5»  /D 

IF(KP(4J.EQ.l)  W9ITE(6,7»  ITER,SRAT,B 
CALL  NULS ( 12, 13,B) 

CALL  SUBVCl 3,9 , iAj 

S*ETA*SRAT 

CALL  A0DS( 14, 20,S) 

CALL  0IVV(2J, 1,201 
S=-3*(CNE+PDIS) 

CALL  PULS( 10, 15,S) 

S=pai S-ONE 

CALL  VULS(20,18,S) 

S=X(3,l)+(ONE-POIS)*SPAT 
CALL  A0DS(18,19,S) 

CALL  SU3V( 19,3,19) 

S=-B*(ONE+POIS) 

CALL  I'ULSC  10,17, S) 

CALL  AOOV( 19, 17,19) 

ITER=ITER+1 

IF(ITSR.GT.KP(3))  GO  TC  200 
CALL  MULS(4,16,B) 

X (16, 1)=0N.E 
CALL  0IVV(19,16,6) 

GO  TO  105 

200  CALL  ^ULS(4 ,7 ,BFAO) 

CALL  MULS(20,8,SR6) 

CALL  MULS(19,9,SRB) 

CALL  ADDV(8,9,9) 

S=SRB/8RA0**2 
CALL  )'ULS(2,10,S) 

CALL  NULS(3 ,11 ,SRB) 

CALL  MULS(1,15,TEEBEE) 

WRITE (6,204 ) 

204  FOPMAT(////) 

WPITE(6,205 ) 

2)5  FORMAT (24X,» RADIUS' ,1 IX, 'THICKNESS ' , 6X, ' GAMM A .OMEGA . SO • , 8X, 
1'  EE.  ALPHA.  TEE'  , BX, ' SI  GMA. '»AOI  AL ',  7X  ,' SIGMA  .C  IRCUMF'  ) 
KSKIP=KP(2) 

00  210  I=1,N,KSKIP 
J=I/KSKIP 

WRITE(6,211  )J,X(7, 1 ), X(15,I ) ,X( 1), I ) ,X(11,I ) ,X( e,I ) ,X(9, 1) 

210  CONTINUE 

211  FORMAT!  no,  1P6E20.5) 

RETURN 

END 


- 12  - 


- 


THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
JBOll  COPY  FURNISHED  TO  DDC  -- 


Listing  of  ancillary  subroutines 


SUBflOLTINe  AD0V(Nl»N2,^3) 

ANCILLARY  SUBROUTINE 

AOC  TWO  VECTCRS  TERM  3Y  TERM 

X(N3,I )=X(N1 »1)+XCN2,I ) 

REALMS  X(23,101  ),S 
COMMON  X,N 
00  1 1=1, N 

1 X(N3, I)=X(Nl,l)+X(N2, I) 

RETURN 

END 

SL/BROUTINE  SUB V (N1 , N2  ,N3 1 
ANCILLARY  SUPROUTINf 
SUBTRACT  TWO  VECTORS  ERM  8Y  TERM 
^ X(N3,I  ) = X(Nl,n-X(M2,n 
REALMS  X(2J,131),S 
COMMON  X,N 
00  1 1=1, N 

1 X(N3, I )=X(N1, I »-X(M2,I ) 

RETURN 


SUBROUTINE  MULV (Nl, N2,N3 ) 

ANCILLARY  SUBROUTINE 

MULTIPLY  TWO  VECTORS  TERM  BY  TERM 

X (N3,I)  = X(N1,I )»X(N2,I  ) 

REAL’^8  X(23,101I,S 
COMMON  X,N 
DC  1 I=1,N 

1 X(N3,I)=X(N1,I)*X(N2,I » 

RETURN 

END 


SUBROUTINE  DIVV (Nl , N2, N3 ) 

ANCILLARY  SUBROUTINE 

DIVIDE  TWO  VECTORS  TERM  BY  TERM 

X(N3,!  ) = X(N1  ,n/X(N2,  ! ) 

REAL*6  X(20,101),S 
COMMON  X,N 

1 x(N3,  i)U(Ni,n /x(N2,n 
RETURN 
END 

SUBROUTINE  ADDS (Nl , N2, S » 

ANCILLARY  SUBROUTINE 

ADC  A SCALAR  TO  EACH  TERM  OF  A VECTOR 
X(N2,n  = X(Nl,IUS 
REALMS  X(20,101),S 
COMMON  X,N 
DO  1 1=1, N _ ^ 

1 X(N2,n=X(Nl,I)+S 
RETURN 
END 
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C 
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SU9IICUTINE  SUBS  (N1.-,N2,') 

ANCILLARY  SUEROUTINP  , 

SUBTRACT  A SCALAR  FROM  EACH  TERM  CF  A VECTOR 
X(N2,I»=X(Nl»I)-S 
REAL*e  X(20,101)»S 
COMMON  X,N 
00  1 I=ltN 
X1N2,I»=X(NI,I)-S 
RETURN 
5NC 


SUBROUTINE  MULS (Nl , N2 , S ) 
ANCILLARY  SUBROUTINE 
MULTIPLY  EACH  TERM  CF 
X(N2,I)=X(Nl,I)*S 
REAL*e  X{20,101)tS 
COMMON  X,N 
DO  L 1=1, N 
X(N2,I)=X(N1,I»*S 
RETURN 
END 

SLBROLTINE  DI VS (Nl , N2 , S ) 
ANCILLARY  SUBROUTINE 
•DIVIDE  EACH  TERM  OF  A 
X(N2,I  ) = X(N1,I)/S 
REAL*8  X(20,101),S 
COMMON  X,N 
00  L 1=1, N 
X(N2,  n = X(Nl,  I ) /S 
RETURN 
END 


A VECTOR  BY  A SCALAR 


VECTOR  BY  A SCALAR 


SUBROUTINE  0UPV(N1,N2) 

ANCILLARY  SUBROUTINE 
DUPLICATE  A VECTOR 
X (N2,I > = X(Nl  ,I > 

REAL*8  X(20,101),S 
COMMON  X,N 
00  I 1=1, N 
1 XlN2,I)=X(Nl,I) 

RETURN 

END 

SUBROUTINE  INTV{Nl,N2) 

ANCILLARY  SUBROUTINE 

INTEGRATE  A VECTOR  USING  MILNE'S  METHOD 
REAL*8  X(20, 101 ),EN,R, add, NINO, NTN0,FIV0,THT0 

COMMON  X,N  i 

EN=N-1 

£N=1.0+0/EN 

R=EN/2.4D+1 

NINa=R*9.D+0 

N'^N0  = R«1.90  + 1 

FIVO=R*5.0+0  ’ 

THT0=R*l.30+l 

X(N2, l)=C.D+0 

X (N2,2)=NINC*X(M1,1 )+NTNO*X(Nl, 2)-F IVC*X(N1,3)+R*X(N1,4) 

NM3=N-3 
00  1 K=1,NM3 

KPl*K-»l  > 

KP2SX  + 2 ) 

KP3  = K-f3  i 

AOD’THTOOI X(N1,KPI)+X(N1,KP2) )-R*( X(NI,K)+X(N1,KP3) I 
1 X(N2,KP2)=X{N2,KP1 )+ACO 

X(N2,N)=X(N2,N-l)+NIN0*X(Nl,N)+NTN0*X(M  ,N-l ) -F  I VC*X  ( Nl  ,N-2  ) J 
l+R*X(M,N-3) 

RETURN 

cNC 
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Appendix  C 


Examples 

The  listing  on  the  next  page  is  of  a MAIN  program  which 
calls  RODISK  twice  to  solve  two  different  problems.  On  am  IBM 
360/67  compile  time  (for  the  MAIN,  RODISK,  and  the  ancillary 
subroutines)  was  12  s,  link  time  was  2 s,  aind  execution  time 
for  both  broblems  was  1.5  s.  In  both  problems  we  used  N = 4l 
and  Iterated  10  times. 

The  first  problem  is  of  type  2 (solid  disk)  with  b = 10, 
t = l/(1.6+.008r^),  V =•  0.3,  “ 120,  o^(b)  = 14000,  and 

BaT  = 25105  + 1300  logg(t)  - r*(233+l6t)  Units  are  inches  and 
pounds.  The  problem  was  made  up  from  the  exact  solution 
a = 9000  + 50r^;  Oq  = + r^(120+l6t) 

r dr 

The  results,  shown  on  paige  I8,  show  evaluations  for  the 

stress  components  which  are  correct  within  0.03  psi  even  though 

convergence  was  complete  to  only  about  4 digits  as  indicated  by 

the  sequence  of  values  above  the  final  tabulation;  these  are 

values  of  a (0)/o„(b)  emd  of  B. 

r r 

The  second  problem  is  of  type  1 with  a » 2.57,  b * 5.15, 
Oj,(a)  - 18205,  0^,(8)  - 22000,  t - 0.1493  T = 60  - 1.6rS 
V » 0.3,  and  Ba  ■ 194.3.  The  units  are  inches  and  pounds. 

Since  the  thickness  variation  is  a power  law,  the  theoretical 
solution  given  in  the  body  hereof  may  be  used  to  obtain  the 
theoretically  exact  solution 

a • -113. 95r*  + 15832rP  - 11708r‘’ 
r 

Oq  ■ 122. 80r*  + 13801rP  + l4310r^ 
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J.  =.  BP.OCK  1 may  1976 

THIS  IS  A MAIN  POOGkAM  TO  TEST  MY  SUBROUTINE  RODISK 
IMPLICIT  REAL#8(A-H,0-Z» 

REALMS  X(20tl31) 

INTEGER  KP(4) 

COMMON  X,N 

COMMON  /ClNE/ARAO,BRAO  ,SRA  ,SR6,TEEB£E,PJIS,KP 
THE  FIRST  PROBLEM  IS  FOP.  A SOLID  DISK 
KP( 1) =2 
KP(2)=4 
KPO)  =10 
KP(4) =1 
BRA0=1.D+1 
SR3=1.40+4 
N = 41 
ENM=N-1 
00  20  1=1, N 
EIM=I-1 

R=EIM*BRAD/ENM 

THICK=l.D+0/( 1 .6D+0+8.D-3^P^*2» 

X(2J,I )=THiCK 

X(2,I  )=120.D+0=BRAO’!'=!‘2/SR3 

W=2. 5 1050+4+1. 30+ B^OLOG { thick )-(R#=2)* (2. 33D+2+1.6D+l#THICK) 
X (3,1 )=W/SR8 
20  CONTINUE 
S=X(2C,NI 
TEE9EE=S 

CALL  CIVS(20,1,S» 

CALL  RODISK 

THE  SECOND  PROBLEM  IS  FOR  AN  ANNULAR  DISK 
TEEBE ==1.4930-1*5.150+0** (-4.20-1) 

ENM=N-1 
KP(1) =1 
ARA0=2. 570+0 
3RA0=5. 150+0 


SRA=1. 82050+4 
SRB=2.20+4 
BETA=5. 024734710-1 
00  30  1=1  ,N 
EIM=I-1 
Y = EIM/£N.M 
X(5,I )=Y 
X(6, I )=Y 

R=ARAC+(BRAD-ARAO)*Y 
X(4,I )=R/eRAO 

T H luK = 1. 4930-1* R**( -4. 20-1) 

X(  1,1  )=''HICK/TEEBEE 
X(2, I ) = BETA 

X (3, 1 ) = (6.0  + l-l .6C+0*R**2  )* 1.9430  + 2/ SRB 
30  CONTINUE 

call  rooisk 
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where  p » .29171  and  q » -1.87171.  The  computer  results  acgree 
with  five  digit  accuracy  even  though  the  sequence  of  iterants 
(A,B)  shows  only  four  digit  convergence. 

The  same  two  problems  were  also  worked  with  N = 101  and  1^ 
Iterations.  The  execution  time  went  from  1.5  s to  3.0  s.  The 
maximum  change  in  any  stress  value  was  0.3  psi.  Prom  these  and 
other  problems  it  may  be  concluded  that  there  are  no  difficulties 
of  accuracy,  computer  storage,  or  execution  time. 
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